This program examines block groups that only partially intersect with community districts.

Load relevant files

# Important abbreviations:
#    bg = block group
#    bl = block
#    cd = community district
#    co = county
#    nbd=neighborhood
#    tr = tract

library(glue)
## Warning: package 'glue' was built under R version 4.1.3
library(magrittr)
## Warning: package 'magrittr' was built under R version 4.1.3
library(sf)
## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.7     v dplyr   1.0.9
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.0.1     v forcats 0.5.1
## Warning: package 'tibble' was built under R version 4.1.3
## Warning: package 'dplyr' was built under R version 4.1.3
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::extract()   masks magrittr::extract()
## x dplyr::filter()    masks stats::filter()
## x dplyr::lag()       masks stats::lag()
## x purrr::set_names() masks magrittr::set_names()
library(dplyr)
path_name <- "../data/"

Contents of bg

load(glue("{path_name}bg.RData"))
glimpse(bg)
## Rows: 1,684
## Columns: 4
## $ bg_id    <chr> "201039819003", "201030714004", "201030711041", "201030712051~
## $ bg_name  <chr> "Block Group 3", "Block Group 4", "Block Group 1", "Block Gro~
## $ bg_area  <dbl> 984975.5, 36602869.0, 18313082.5, 26007918.3, 11448660.8, 243~
## $ geometry <POLYGON [°]> POLYGON ((-94.94756 39.3373..., POLYGON ((-95.18756 3~

Contents of bl

load(glue("{path_name}bl.RData"))
glimpse(bl)
## Rows: 39,903
## Columns: 4
## $ bl_id    <chr> "200910518012010", "200910514003006", "200910515002005", "200~
## $ bl_name  <chr> "Block 2010", "Block 3006", "Block 2005", "Block 1009", "Bloc~
## $ bl_area  <dbl> 38500.823, 51711.677, 104537.047, 73933.762, 95854.528, 19934~
## $ geometry <POLYGON [°]> POLYGON ((-94.64213 38.9774..., POLYGON ((-94.63966 3~

Contents of cd

load(glue("{path_name}cd.RData"))
glimpse(cd)
## Rows: 59
## Columns: 4
## $ cd_id    <dbl> 106, 108, 113, 102, 129, 116, 114, 101, 105, 103, 107, 109, 1~
## $ cd_name  <chr> "East Side", "Old Northeast", "Greater Downtown", "Blue Valle~
## $ cd_area  <dbl> 281713850, 118237219, 181920811, 216842678, 412671242, 264282~
## $ geometry <POLYGON [°]> POLYGON ((-94.52337 39.0941..., POLYGON ((-94.50777 3~

Contents of co

load(glue("{path_name}co.RData"))
glimpse(co)
## Rows: 220
## Columns: 19
## $ STATEFP  <chr> "29", "29", "20", "20", "20", "29", "29", "20", "20", "29", "~
## $ COUNTYFP <chr> "083", "011", "073", "043", "157", "103", "117", "039", "147"~
## $ COUNTYNS <chr> "00758496", "00758460", "00485003", "00484991", "00485042", "~
## $ GEOID    <chr> "29083", "29011", "20073", "20043", "20157", "29103", "29117"~
## $ NAME     <chr> "Henry", "Barton", "Greenwood", "Doniphan", "Republic", "Knox~
## $ NAMELSAD <chr> "Henry County", "Barton County", "Greenwood County", "Donipha~
## $ LSAD     <chr> "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "~
## $ CLASSFP  <chr> "H1", "H1", "H1", "H1", "H1", "H1", "H1", "H1", "H1", "H1", "~
## $ MTFCC    <chr> "G4020", "G4020", "G4020", "G4020", "G4020", "G4020", "G4020"~
## $ CSAFP    <chr> NA, NA, NA, "312", NA, NA, NA, NA, NA, "312", NA, NA, NA, NA,~
## $ CBSAFP   <chr> NA, NA, NA, "41140", NA, NA, NA, NA, NA, "47660", NA, "21380"~
## $ METDIVFP <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N~
## $ FUNCSTAT <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "~
## $ ALAND    <dbl> 1805054746, 1533351028, 2961143636, 1019105691, 1857998463, 1~
## $ AWATER   <dbl> 91657543, 12152201, 24145021, 12424174, 7951912, 7307094, 161~
## $ INTPTLAT <chr> "+38.3864909", "+37.5007989", "+37.8793472", "+39.7885021", "~
## $ INTPTLON <chr> "-093.7926278", "-094.3440893", "-096.2417321", "-095.1472253~
## $ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-93.51565 3..., MULTIPOLYGON (((~
## $ label    <chr> "Henry", "Barton", "Greenwood", "Doniphan", "Republic", "Knox~

Contents of cd-intersections

load(glue("{path_name}cd-intersections.RData"))
glimpse(bg_cd_intersection)
## Rows: 1,236
## Columns: 3
## $ bg_id      <chr> "290950161001", "290950166002", "290950060001", "2909500340~
## $ cd_id      <dbl> 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106,~
## $ bg_prop_in <dbl> 0.8919314491, 1.0000000000, 1.0000000000, 1.0000000000, 0.3~
glimpse(bl_cd_intersection)
## Rows: 15,748
## Columns: 5
## $ tr_id      <chr> "29095002300", "29095006300", "29095002200", "29095002200",~
## $ bg_id      <chr> "290950023002", "290950063002", "290950022003", "2909500220~
## $ bl_id      <chr> "290950023002014", "290950063002017", "290950022003001", "2~
## $ cd_id      <dbl> 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106,~
## $ bl_prop_in <dbl> 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 0.997, 0.001, 0.0~

Contents of cd-weights

load(glue("{path_name}cd-weights.RData"))
glimpse(bg_counts)
## Rows: 1,236
## Columns: 6
## Groups: bg_id [776]
## $ bg_id     <chr> "200910500001", "200910500002", "200910500003", "20091050100~
## $ cd_id     <dbl> 209, 209, 209, 209, 211, 209, 211, 211, 112, 112, 112, 111, ~
## $ people_in <dbl> 1, 2, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, ~
## $ units_in  <dbl> 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ~
## $ people    <dbl> 744, 695, 751, 1313, 1313, 1196, 643, 630, 778, 961, 1210, 9~
## $ units     <dbl> 341, 311, 377, 648, 648, 549, 315, 444, 287, 335, 461, 430, ~
glimpse(bg_list)
##  chr [1:776] "290950161001" "290950166002" "290950060001" "290950034004" ...
glimpse(bl_counts)
## Rows: 41,220
## Columns: 6
## $ bg_id     <chr> "290950023002", "290950063002", "290950022003", "29095002200~
## $ bl_id     <chr> "290950023002014", "290950063002017", "290950022003001", "29~
## $ cd_id     <dbl> 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, 106, ~
## $ prop_in   <dbl> 1.000000000, 1.000000000, 1.000000000, 1.000000000, 1.000000~
## $ people_in <dbl> 62, 0, 22, 26, 42, 49, 43, 0, 0, 20, 5, 56, 58, 48, 47, 13, ~
## $ units_in  <dbl> 18, 0, 9, 18, 22, 22, 14, 0, 0, 10, 4, 12, 19, 24, 28, 4, 48~
glimpse(bl_list)
## Rows: 17,547
## Columns: 2
## $ bg_id <chr> "200910533013", "200910533024", "200910533024", "200910534291", ~
## $ bl_id <chr> "200910533013008", "200910533024013", "200910533024008", "200910~

Subset co to seven counties

clist <- c(
  "20091",
  "20103",
  "20209",
  "29037",
  "29047",
  "29095",
  "29165")
co                                     %>%
  filter(GEOID %in% clist)             -> co
glimpse(co)
## Rows: 7
## Columns: 19
## $ STATEFP  <chr> "29", "29", "20", "20", "29", "29", "20"
## $ COUNTYFP <chr> "037", "165", "091", "209", "095", "047", "103"
## $ COUNTYNS <chr> "00758473", "00758537", "00485010", "00485065", "00758502", "~
## $ GEOID    <chr> "29037", "29165", "20091", "20209", "29095", "29047", "20103"
## $ NAME     <chr> "Cass", "Platte", "Johnson", "Wyandotte", "Jackson", "Clay", ~
## $ NAMELSAD <chr> "Cass County", "Platte County", "Johnson County", "Wyandotte ~
## $ LSAD     <chr> "06", "06", "06", "06", "06", "06", "06"
## $ CLASSFP  <chr> "H1", "H1", "H1", "H6", "H1", "H1", "H1"
## $ MTFCC    <chr> "G4020", "G4020", "G4020", "G4020", "G4020", "G4020", "G4020"
## $ CSAFP    <chr> "312", "312", "312", "312", "312", "312", "312"
## $ CBSAFP   <chr> "28140", "28140", "28140", "28140", "28140", "28140", "28140"
## $ METDIVFP <chr> NA, NA, NA, NA, NA, NA, NA
## $ FUNCSTAT <chr> "A", "A", "A", "C", "A", "A", "A"
## $ ALAND    <dbl> 1804223313, 1087283058, 1226679688, 392739381, 1565698757, 10~
## $ AWATER   <dbl> 14963188, 16949451, 16319024, 11801597, 30621016, 28508804, 1~
## $ INTPTLAT <chr> "+38.6464737", "+39.3786900", "+38.8839065", "+39.1153842", "~
## $ INTPTLON <chr> "-094.3545467", "-094.7614765", "-094.8223295", "-094.7630866~
## $ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-94.11966 3..., MULTIPOLYGON (((-94.65596 3..~
## $ label    <chr> "Cass", "Platte", "Johnson", "Wyandotte", "Jackson", "Clay", ~

Plot functions

# These functions plot various census, 
# community, or neighborhood  geographies.

add_red_border <- function(mapx, geo) {
  mapx                                    +
    geom_sf(
      data=geo,
      aes(),
      fill=NA,
      color="darkred")
}

add_green_interior <- function(mapx, geo) {
  mapx                                    +
    geom_sf(
      data=geo,
      aes(),
      fill="lightgreen",
      color=NA)
}

Get id values for all community districts

cd                                     %>%
  tibble                               %>%
  distinct(cd_id)                      %>%
  arrange(cd_id)                       %>%
  pull(cd_id)                           -> cd_list

Draw an overview map

co_points <- st_point_on_surface(co)
## Warning in st_point_on_surface.sf(co): st_point_on_surface assumes attributes
## are constant over geometries of x
## Warning in st_point_on_surface.sfc(st_geometry(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
co_coords <- as.data.frame(st_coordinates(co_points))
co_coords$NAME <- co$NAME

ggplot()+
  geom_sf(
    data=co, 
    color="darkred", 
    fill=NA)                            +
  geom_text(
    data=co_coords, 
    aes(X, Y, label=NAME))              +
  xlab("Longitude")                     +
  ylab("Latitude")                      -> overview_map
plot(overview_map)

High level view

for (i_cd in cd_list) {
  # Draw community district in green atop
  # counties outlined in red
  cd_subset <- filter(cd, cd_id==i_cd) 
  ti <- cd_subset$cd_name
  overview_map                          +
    ggtitle(ti)                         +
    geom_sf(
      data=cd_subset,
      aes(),
      fill="lightgreen",
      color=NA)                         -> map1
  plot(map1)

  # draw community districts in green with
  # block groups outlined in red.
  bg_cd_intersection                   %>%
    filter(cd_id==i_cd)                %>%
    filter(bg_prop_in > 0.10)          %>%
    pull(bg_id)                         -> bg_list
  bg                                   %>%
    filter(bg_id %in% bg_list)          -> bg_subset

  ggplot(cd_subset)                     +
    ggtitle(ti)                         +
    xlab("Longitude")                   +
    ylab("Latitude")                    +
      geom_sf(
        data=cd_subset,
        aes(),
        fill="lightgreen",
        color=NA)                       +
    geom_sf(
      data=bg_subset,
      aes(),
      fill=NA,
      color="darkred")                  -> map2
  plot(map2)

  # draw block groups and blocks outlined
  # in red atop community district in green.
  for (i_bg in bg_list) {
    bg_cd_intersection                 %>%
      filter(cd_id==i_cd)              %>%
      filter(bg_id==i_bg)              %>%
      pull(bg_prop_in)                  -> p
    if (p < 0.1) next
    if (p > 0.9) next
    bl_subset <- filter(bl, str_sub(bl_id, 1, 12)==i_bg)
    ggplot(cd_subset)                   +
      ggtitle(ti)                       +
      xlab("Longitude")                 +
      ylab("Latitude")                  +
    geom_sf(
      data=cd_subset,
      aes(),
      fill="lightgreen",
      color=NA)                         +
    geom_sf(
      data=bl_subset,
      aes(),
      fill=NA,
      color="darkred")                  -> map3
    plot(map3)
  }
}

for (i_cd in cd_list[5:7]) {
  next
  bg_n1 <- dim(bg_subset)[1]
  bg_n2 <- sum(bg_subset$bg_prop_in==1)
  bg_n3 <- sum(bg_subset$bg_prop_in > hi & bg_subset$bg_prop_in != 1)
  bg_n4 <- sum(bg_subset$bg_prop_in > lo & bg_subset$bg_prop_in <= hi)
  bg_n5 <- sum(bg_subset$bg_prop_in <= lo)
  
  bg_message <- glue(
    "{bg_n1} block groups touch or intersect community district\n",
    "{bg_n2} block groups lie entirely inside community district\n",
    "{bg_n3} block groups lie mostly inside community district\n",
    "{bg_n4} block groups lie partially inside community district\n",
    "{bg_n5} block groups lie mostly outside community district\n")
  cat(bg_message)
  plot_cd(cd, i_cd, bg, bg_subset$bg_id, "block groups touching or intersecting")
  plot_cd(cd, i_cd, bg, bg_subset$bg_id, "block groups touching or intersecting")
    for (i_bg in bg_subset$bg_id) {
    bg_cd_intersection                   %>%
      filter(cd_id==i_cd)                %>%
      filter(bg_id==i_bg)                %>%
      pull(bg_prop_in)                    -> bg_prop
    if (bg_prop < lo) next
    if (bg_prop > hi) next
    plot_cd(cd, i_cd, bg, i_bg)
    plot_bg(cd, i_cd, bg, i_bg, bl)
    bl_cd_intersection                   %>%
      filter(cd_id==i_cd)                %>%
      filter(bg_id==i_bg)                %>%
      select(bl_id, bl_prop_in)           -> bl_subset
    bl_n1 <- dim(filter(bl_cd_intersection, bg_id==i_bg))[1]
    bl_n2 <- sum(bg_subset$bg_prop_in==1)
    bl_n3 <- sum(bg_subset$bg_prop_in > hi & bg_subset$bg_prop_in != 1)
    bl_n4 <- sum(bg_subset$bg_prop_in > lo & bg_subset$bg_prop_in <= hi)
    bl_n5 <- sum(bg_subset$bg_prop_in <= lo)
    bl_n6 <- bl_n1 - bl_n2 - bl_n3 - bl_n4 - bl_n5
    bl_message <- glue(
      "{bl_n1} blocks in the block group\n",
      "{bl_n2} blocks lie entirely inside community district\n",
      "{bl_n3} blocks lie mostly inside community district\n",
      "{bl_n4} blocks lie partially inside community district\n",
      "{bl_n5} blocks lie mostly outside community district\n",
      "{bl_n6} blocks lie entirely outside community district\n")
    cat(bl_message)
    for (i_bl in bl_subset$bl_id) {
      bl_cd_intersection                 %>%
        filter(cd_id==i_cd)                %>%
        filter(bg_id==i_bg)                %>%
        filter(bl_id==i_bl)                %>%
        pull(bl_prop_in)                    -> bl_prop
      if (bl_prop > hi) next
      if (bl_prop < lo) next
      plot_bl(cd, i_cd, bg, i_bg, bl, i_bl)
    }
    bg_counts                            %>%
      filter(cd_id==i_cd)                %>%
      filter(bg_id==i_bg)                %>%
      data.frame                         %>%
      print
    bl_counts                            %>%
      filter(cd_id==i_cd)                %>%
      filter(bg_id==i_bg)                %>%
      select(-bg_id)                     %>%
      data.frame                         %>%
      print
  }
}